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Abstract 

We study fluid-fluid phase separation in a colloid-polymer mixture adsorbed in a colloidal porous 
matrix close to the 9 point. For this purpose we consider the Asakura-Oosawa model in the presence 
of a quenched matrix of colloidal hard spheres. Wc study the dependence of the demixing curve 
on the parameters that characterize the quenched matrix, fixing the polymer-to-colloid size ratio 
to 0.8. We find that, to a large extent, demixing curves depend only on a single parameter /, 
which represents the volume fraction which is unavailable to the colloids. We perform Monte Carlo 
simulations for volume fractions / equal to 40% and 70%, finding that the binodal curves in the 
polymer and colloid packing-fraction plane have a small dependence on disorder. The critical point 
instead changes significantly: for instance, the colloid packing fraction at criticality increases with 
increasing /. Finally, we observe for some values of the parameters capillary condensation of the 
colloids: a bulk colloid-poor phase is in chemical equilibrium with a colloid-rich phase in the matrix. 



PACS: 61.25.Hq, 82.35.Lr 



I. INTRODUCTION 



The study of the fluid phases in mixtures of colloids and nonadsorbing neutral poly- 
mers has become increasingly important in recent years; see Refs. [l|-6| for recent reviews, 



Refs 



.0 



16| for experiments, and Refs. 17H4d| for theoretical investigations. These sys- 



tems show a very interesting phenomenology, which only depends to a large extent on the 
nature of the solvent and on the ratio q = Rg/Rc, where Rg is the radius of gyration of 
the polymer and Rc is the radius of the colloid. Experiments and numerical simulations 
indicate that polymer-colloid mixtures have a solid colloidal phase for large enough colloidal 
concentrations and a corresponding fluid-solid coexistence. Much less obvious is the pres- 
ence of a fluid-fluid coexistence of a colloid-rich, polymer-poor phase (colloid liquid) with a 
colloid-poor, polymer-rich phase (colloid gas). Extensive theoretical and experimental work 
indicates that such a transition occurs only if the size of the polymers is sufficiently large. 



i.e. for q > q*, where 



38| g* ^ 0.3-0.4. 



At least qualitatively, many aspects of the behavior of colloid-polymer suspensions can be 



41 



42l | , which gives a coarse-grained 



understood by using the Asakura-Oosawa (AO) model 
description of the mixture. The polymers are treated as an ideal gas of point particles of 
radius Rp (which is usually identified with the radius of gyration) which interact with the 
colloids by means of a simple hard-core potential. This model is extremely crude since it 
ignores the polymeric structure and polymer-polymer repulsion, which is relevant in the 
good-solvent regime. Nonetheless, it correctly predicts polymer-colloid demixing as a result 
of the entropy-driven effective attraction (depletion interaction) between colloidal pairs due 



to the presence of the polymers 



17H19, 



21 



23|-l25, 



31 



32l |. It is not, however, quantitatively 



predictive for polymers in the good-solvent regime. For instance, at a given colloid packing 
fraction, the AO model predicts the binodal curve to be at a polymer volume fraction which is 
significantly lower than that observed experimentally. In order to reproduce the experimental 
results one can use realistic atomistic models for the polymers, but this is a very difficult 
task from a numerical point of view. In the colloid regime g < 1, it is much easier, and still 



provides good results, to use coarse-grained models in which polymers are mode 
particles (as in the AO model) interacting with repulsive soft pair potentials 24 



ed as point 



26 



33 



4Q|, 



which have either a phenomenological origin or are derived by means of exact coarse-graining 
procedures. Nonetheless, at least for g < 1 (colloidal regime), the AO model is expected 
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to provide quantitatively correct results for colloid-polymer solutions close to the 6 point. 
Indeed, in this regime polymers show an approximate ideal behavior and can be described 
quite reasonably as noninteracting random walks, as does the AO model 43[ |. 

In this paper we wish to study the demixing of colloid-polymer mixtures in porous ma- 
terials, which are characterized by a highly interconnected porous structure. They have 
important technological applications, for instance in catalysis and gas separation and pu- 



rification 



44j |. Examples are the Vycor glasses, in which pore sizes range from 1 nm to 100 



nm, and high-porosity systems like silica gels (xerogels and aerogels), which are produced 



by means of silica sol- 
been studied in Refs. 



ge 



processes. AO colloid-polymer mixtures in a porous matrix have 
45|-|49| by means of density-functional theory, integral equations, and 



Monte Carlo (MC) simulations. The nature of the critical transition has been fully clarified 



46N49l|: if obstacles are random and there is a preferred affinity of the quenched obstacles 



to one of the phases, the transition is in the same universality class as that occurring in the 
random- field Ising model, in agreement with a general argument by de Gennes 50|. If these 
conditions are not satisfied, standard Ising or randomly dilute Ising behavior is observed 



instead, see Refs. 



51 



, . On the other hand, little is known on how demixing is infiuenced 
by the amount of disorder and by its nature (for a polymer matrix some results for the 
critical-point behavior as a function of the amount of disorder are reported in Ref. 46|). In 
this paper porosity is introduced by considering a quenched matrix of hard spheres of radius 
-Rdis- We will compute the binodal curves in terms of the polymer and colloid packing frac- 
tions for different ratios Rdis/Rc and for different disorder concentrations with the purpose 
of determining how these parameters affect the location of the demixing transition and of 
the critical (second-order) transition point. We will not instead perform a detailed study 
of the q dependence and we shall set q = 0.8 as in Ref. [481]. This work complements the 
results of Ref. |45|, which instead studied the q dependence for a single value of Rdis/Rc, 
Rdis/Rc = 1, and of the disorder concentration. 

The paper is organized as follows. In Sec. [Ill we discuss the model and the relevant 
variables. In Sec. IIIII we present our numerical results. Our conclusions are presented in 
Sec. IIVI In App. |A]we present some details on the MC calculation. 
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II. THE MODEL 



In the AO model polymers and colloids are modelled as spheres of radii Rp and Rc, 
respectively. We assume hard-sphere interactions between colloid and colloid-polymer pairs; 
the pair potentials are given by 

for r < 2i?c, 
for r > 2Rc, 

for r < Rc + Rp, 
for r > Rc + Rp, 

(1) 

where r is the center-to-center distance. We consider a cubic box of size L and we characterize 
the thermodynamic phases in terms of the packing fractions 

Vp- , l^j 

where Np and Nc indicate the number of polymers and of colloids in the box, respectively. 

The phase behavior of the AO model has been extensively studied. It strongly depends 
on the polymer-to-colloid size ratio q = Rp/ Rc- For small values of q the demixing transition 
is unstable and only the fluid-solid transition occurs. Fluid- fluid demixing occurs jl|, H, 38 1 




for q > 0.3-0.4. In this work we have not investigated the q dependence of the binodal curve, 
since our main objective is the analysis of the role of quenched disorder. We have thus fixed 



q = 0.8, as in Ref. [48||, at the boundary between the colloid and the protein regimes. 

Disorder has been introduced by considering a colloidal quenched matrix which has a 
hard-sphere interaction both with the colloids and the polymers. In practice, we choose a 
disorder concentration Cdis and randomly distribute A^^dis = CdisL^ nonoverlapping spheres of 
radius Rdis in the box. The position of these spheres is assumed to be fixed (quenched). 
Colloids and polymers can only move outside the quenched matrix, which means that the 
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spheres belonging to the matrix and the freely moving particles interact with pair potentials 

Mc,dis(^) = 



CO for r < Rc + Rdis, 

for r > + Rdis, 



oo for r < _Rp + _Rdis, 

Up^isir) = { (3) 

ior r > Rp + i?dis- 



Note that the matrices considered here are different from those discussed in Refs. |46|, |4j 
The main difference is that here the matrix consists in hard spheres that cannot intersect 
each othe. ,we .a,„e it coUoida. ™atn.). 0„ the othe. ha.d. Ref. Q the ,„at.x 
spheres are soft and can freely overlap, as if they were an ideal gas (hence the name polymer 
matrix). Second, in those works, for a given choice of Cdis, the number A^dis is not fixed, 
but obtained from a Poissonian distribution with mean value CdisL^- This second difference 
should not be important in the infinite- volume limit, since it entails density fluctuations of 
order l/L^^^, which vanish as L — cxd. 

In the simple model we consider, disorder is characterized by two parameters, the reduced 
concentration c = Cdis-Rj? ratio Rdis/ Rc- However, c does not directly characterize the 

free space available to the colloids and to the polymers. We shall use instead the effective 
matrix filled-space ratio /, which is defined as follows. Consider the region TZ in which the 
(centers of the) colloids are allowed: 

n = {v:\v-Vi\>R, + Rdis, for all 1 < i < Ndis }, (4) 

where is the position of the i-th hard sphere belonging to the matrix. If Vji is the volume 
of the region TZ, we define 

/ = 1 - -jr, (5) 

where [Vji] is the average of Vji over the different matrix realizations. Note that, for large 
values of L, [Vti] is essentially independent of the matrix realization, a property known as self- 
averaging. The parameter / represents the volume fraction that is unavailable to the colloids 
due to the presence of the random matrix and can easily be determined by computing the 
probability of inserting a colloid in the otherwise empty matrix. In a completely analogous 
way we can define /poi, which characterizes the volume fraction unavailable to polymers. Of 
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FIG. 1: Two-dimensional systems with L/Rc = 50 and / = 0.5. On the left we take Rdis/Rc = 0.1, 
while on the right Rdis/Rc = 3.0. The disorder packing fractions TrciJ^jg are 0.0057 and 0.292 in 
the two cases, respectively. The gray circles of radius Rc + Rdis correspond to the colloid-excluded 
region (depletion region) around each sphere of the quenched matrix: the centers of the colloids 
can only belong to the white region. We also draw a single colloid (black) to show the length scale. 

course, / > /poi in the colloid regime in which g < 1, while / < /poi in the opposite, protein 
regime. 

It is interesting to understand qualitatively how the disorder distribution changes with 
i?dis at fixed /. In Fig. [T]we show the matrix for / = 0.5 and two values of i^dis, -Rdis/ -Rc = 0.1 
and Rdis/Rc = 3. To make the figure more clear, we consider a two-dimensional system, that 
is a matrix of nonoverlapping disks on a square of area L^. It is evident that the topology 
of the matrix is quite different. For large Rdis/Rc the free volume available to the colloids 
consists in large empty regions connected by narrow channels. This is the case of a porous 
material with big interconnected pores. On the other hand, for Rdis/Rc small, pores are 
significantly smaller and the topology of the network is more complex. 

In order to have demixing, the parameter / cannot be arbitrarily close to 1, but should 
satisfy / < /*, where /* is related to the percolation threshold of the region TZ in which 
colloids can move. For f > f* the space TZ divides in disconnected finite regions and thus 
no phase transition is possible. The exact value of /* is unknown. However, the arguments 
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TABLE I: Estimates of the reduced concentration c = Cdis-Rc ^'^d of the disorder packing fraction 
^dis = 47ri?jjgCdis/3 for two values of / and several values of Rdis/Rc- 

f = 0.40 / = 0.70 

Rdis/Rc c r]dis c rjdi 



IS 



0.005 0.120 6.3 • 10-8 0.283 1.5 • lO^^^ 

0.01 0.118 4.9-10"'^ 0.279 1.2-10"^ 

0.02 0.115 3.9-10-6 0.271 9.1 • 10"^ 

0.05 0.105 5.5 • 10-5 0.248 1.3 • 10-"^ 

0.1 0.0915 3.8-10-^ 0.215 9.0-10-^ 

0.2 0.0700 2.3 • 10-3 0.164 5.5 • 10-^ 

0.4 0.0431 0.0116 0.0972 0.026 

0.6 0.0280 0.0254 0.0609 0.055 

1.0 0.0136 0.057 0.0278 0.116 



of Ref. |53| suggest 



/* ^ 0.85. (6) 

For the same reasons — polymers should be able to move in the whole space — the polymer 
parameter /poi must satisfy /poi < /* in order to observe coexistence. 

In this paper we shall perform simulations for two values of /, / = 0.40 and / = 0.70, 
the latter being quite close to the threshold /*, and for q = 0.8, so that /poi < /. In 
Table [T] we report the reduced concentration c and the disorder packing fraction r^dis = 
47r_R^;gC/3 for several values of Rdis/Rc- First, we observe that c converges to a finite positive 
constant as Rdis/Rc 0. This result is quite easy to understand. If Rdis/Rc ^ 1, the pair 
potentials ([3]) become essentially independent of i?dis- Hence, the density becomes essentially 
independent of i?dis for -Rdis small. In the opposite limit Rdis/Rc ^ 1, the potentials become 
essentially independent of Rc. Hence, in this limit / converges to the disorder packing 
fraction ?7dis = 4:TTR^^^Cdis/3- For instance, for ?7dis = 0.30, we obtain / = 0.51,0.40,0.32 for 



Rdis/Rc = 5,10,50. Since a hquid hard-sphere phase exists only up [54] to r/ ~ 0.49, for 
large Rdis/Rc, the matrix may belong to different hard-sphere phases, while still satisfying 
the condition f < f*. 
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III. RESULTS 



A. Monte Carlo simulation 

In this work we investigate the effect of disorder on the fluid-fluid binodals for q = 0.8. 
We perform simulations in the absence of the porous matrix — our results are consistent 
with those of Refs. 32|, |33| — and for two values of /, / = 0.4 and / = 0.7 [note that 
c(/ = 0.7) ~ 2c(/ = 0.4)], in cubic boxes with L/Rc = 16 and 20. In order to obtain 
quenched averages we consider 200-400 matrix realizations for each / and -Rdis- 

For each value of / we consider a few values of Rdis/Rc- For f = 0.4 we present results 
for Rdis/Rc = 0.2,0.6,1.0 (/poi = 0.26, 0.29, and 0.31, respectively), while for / = 0.7, we 
use Rdis/Rc = 0.2 and 1.0 (/poi = 0.50 and 0.57, respectively). It should be noted that we 
are limited by computer power in further decreasing or increasing the ratio Rdis/Rc- Indeed, 
if we further decrease the ratio, the disorder density increases, see Table [U and so does 
the number of matrix particles and the computational work. On other hand, if we increase 
Rdis/Rc beyond 1, we should consider quite large systems in order to avoid large size effects, 
which is unfeasible with our present computer power. 

In order to determine the coexistence curves we perform a grand-canonical simulation. 
The grand partition sum for each disorder realization is 

EiV, z„ z,) = J2 ^p'^c'QiV, N„ iV,), (7) 

where Q{V, Np, Nc) is the configurational partition function of a system of Np polymers and 
Nc colloids in a volume V, and Zp and Zc are the corresponding fugacities. In Eq. ([71) we 
normalize Q{V, Np, Nc) so that Q{V, 1, 0) = Q{V, 0, 1) = V, hence Zp and Zc are dimensionful 
parameters. We quote our results in terms of the dimensionless combinations z^R^ and 

v; = -j^pRl- (8) 

The quantity Vjp represents the polymer reservoir packing fraction. 

In the presence of a first-order transition, standard local algorithms are unable to sample 
correctly both phases in the simulation. We therefore combine the grand-canonical algorithm 



with the umbrella sampling and the simulated-tempering method [55|, |56|], as discussed in 



App. |Xl Insertions and deletions of colloids and polymers are performed by using the cluster 



moves introduced by Vink and Horbach 



31 



57|. 



B. Quenched coexistence curve 



The main purpose of this work is the determination of the disorder-averaged coexistence 
curve. In order to define it precisely, let us define the disorder- averaged colloid and polymer 
numbers 

d 

Nc,ay{V, Zp, Zc) = z^^ [lnE;(V, Zp, Zc)] , 

OZc 

d 

A^p,av( V, Zp, Zc) = Zp— [In S( V, Zp, Zc)] , (9) 

azp 

where [•] indicates the average over the matrix realizations. In the presence of first-order 
transitions, there is a line Zc — z*{zp) in the {zp, Zc) plane where these two functions become 
discontinuous in the infinite- volume limit. In other words, for Zp > Zp^crit we have 

lim lim A^c,av(V, Zp, z*{zp) + e)/V = Cc,iiq, 

e->0+ V—^oo 

hm lim A^p,av(V", Zp, z*(zp) + e)/V = Cp^uq, 

lim lim iVc,av(V, Zp, Z*{Zp) - €)/V = Cc,gas, 
e^0+ y — >oo 

lim lim iVp,av(V, Zp, z*{zp) - e)/V = Cp,gas. (10) 

e— >-0+ V-^oo 

The pair Cp^uq, Cc,iiq gives the polymer and colloid concentrations in the colloid-liquid phase 
at coexistence, while Cp gas and Cc,gas correspond to the colloid-gas polymer-rich phase. 

In the MC simulations the position of the demixing curve can be determined by studying 
the disorder averaged histograms of Nc and Np, which are defined as 

hc,3,ve{Nc,0, Zp, Zc) = {HNc,Nc,o))GC,Zr„z. ' (^1) 

hp,s.ye{Np,o, Zp, Zc) = [{SiNp, Np,o))gc,z„z] ' (12) 

where 5{x,y) is the Kronecker's delta [5{x,x) — 1, S{x,y) — ior x ^ y] and {■)gczpZc 
is the grand-canonical ensemble average. In the two-phase region the histograms show a 
double-peak structure. In order to obtain z* at fixed Zp in a finite volume, several different 
methods can be used. We followed two different recipes, the equal-area and the equal-height 
methods. In the first case, we define z* as the value of the colloid fugacity at which the area 
below the two peaks is equal. For instance, if we consider the colloid-number distribution, 
we first compute the position A^min of the minimum between the two peaks and then require 
z* to be the value of the colloid fugacity at which 

^ ^ ^c,ave(-^C) -^p) ) ^ ^ ^c,ave(-^C) ^p; ) • i^^) 
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Equivalently, one can use the polymer distribution hp^^veiNc, Zp, Zc). In the second method 
we identify z* as the value of the fugacity at which the two peaks have the same height. 
Once z* has been obtained, the colloid and polymer number at the transition are defined as 
the positions of the maxima of the histograms. 

Since we have two different histograms to analyze, we obtain two different estimates of 
the colloid fugacity at coexistence: an estimate z*{c) is obtained from the analysis of the 
colloid- number histograms, while z* {p) is obtained from the analysis of the polymer- number 
histograms. For R^is/ Rc = 0.2,0.6 the two estimates are quite close and provide consistent 
estimates of the colloid and polymer packing fractions at coexistence, although the equal- 
area method is more thermodynamically consistent. Indeed, the differences \z*{c) — z*{p)\ 
computed with the equal-area method are always smaller than those computed with the 
second prescription. For R^is/Rc = 1-0, we have been unable to apply the area method. The 
difficulties can be understood by looking at Fig. [21 where we report the colloid histograms 
for / = 0.7 and for the largest value of Zp we consider, rjp ^ 1.82 {zpR^ = 0.85). While the 
colloid-liquid peak is quite narrow, the colloid-gas peak is very broad and therefore condition 
( ITSj) is satisfied only when the colloid-gas peak is barely visible. However, in this case the 
definition of A^min is ambiguous and thus z* is determined with large uncertainty. In some 
cases, it is even impossible to satisfy the equal-area condition. Thus, for Rdis/Rc = 1-0 we 
only use the equal-height method. Note that the two methods should give identical results 
in the infinite-volume limit. Hence, the difficulties we observe indicate that for this value 
of the parameters finite-size effects are important. The analyses reported in the following 
sections confirm these findings. It is interesting to note that, at variance with what happens 
in the bulk, in the presence of randomness the order parameter distribution shows two well- 
separated nonoverlapping peaks even at the critical point (see Refs. for a discussion 
in the present context). Thus, it is also possible that the difficulties we observe for some 
values of the parameters are related to the fact that they belong to the one-phase region, 
even if the finite-size colloid and polymer histograms are bimodal. 

To give an idea of the performance of the two methods, we report the results for / = 0.4, 
Rdis/Rc = 0.6, rip ^ 1.24 {zpR^ = 0.58), and L/Rc = 16, see Fig. [31 The equal-area method 
gives 

z^Rl^ 130.2 r/,,gas ~ 0.041 r/,,uq ^ 0.294, (14) 
from the analysis of the colloid distribution. The analysis of the polymer distribution gives 
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FIG. 2: Colloid histogram /ic,ave for LjRc = 16, / = 0.7, i?dis/^c = 1-0, % k, 1.82 {z^Rl = 0.85), 
and some values of Zc close to z*. 

the same estimate of z* (we have data for several values of Zc with step Azc = 0.1). If we 
apply instead the equal-height method we obtain z*{c)Rl ~ 129.3 and z*{p)R^ ^ 127.5 from 
the two distributions. The colloid packing fractions at coexistence are therefore 

?3 



z*(c)R 



129.3 



c,gas 



(15) 



0.039 r/c,iiq ^ 0.294, 
z:(p)i?3 = 127.5 r/c,gas ~ 0.038 r^.^iq ^ 0.292. 

The results are very close with each other and consistent with those reported in Eq. ([1] 
Similar conclusions are obtained for the polymer packing fractions at coexistence. For / = 
0.7, Rdis/Rc = 1-0, V? ~ 1-82 {zpRl = 0.85), and L/R^ = 16, the case reported in Fig. El we 
obtain z*{c)Rl ~ 454 and z*{p)Rl ~ 440 from the two distributions (equal-height method). 
At coexistence we find then 

z:(c)R^.=A5A 



0.096 



VcM 



0.257, 



(16) 



'7c, gas 

z*{p)Rl = UO r/e,gas ~ 0.100 r/e,iiq ~ 0.258. 

Even though the estimates of the coexistence colloid fugacity differ somewhat, the two 
estimates of the colloid packing fractions at coexistence are quite close. 



C. Sample-to-sample fluctuations 



It is interesting to understand how the results obtained from the sample average compare 
with those that would be obtained by determining the coexisting phases for each disorder 
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FIG. 3: Polymer histogram /ip,ave and colloid histogram /ic,avc for L = 16i?c, / = 0.4, Rdis/Rc = 0.6, 
rfp 1.24 {zpR^ = 0.58), and several fugacities ZcR^- The thicker curve corresponds to the 
coexistence fugacity obtained by using the equal-area prescription {z*R^ = 130.2), while the other 
two curves correspond to the estimates z*{c)R^ = 129.3 and z*{p)R^ = 127.5 obtained by using 
the equal-height method. 

realization. In Fig H] we report the distributions of the colloid and polymer packing fractions 
at coexistence computed from each disorder configuration. The distributions of the packing 
fractions corresponding to the colloid-liquid phase are very narrow and are centered at 
the value obtained from the analysis of the average distributions. On the other hand, the 
distributions for the colloid-gas phase are broad, especially for _Rdis/-Rc = 0.6 and 1. Since the 
broadness of the distribution is a finite-size effect — we expect the width of the distributions 
to scale as 1/ a/L as L — oo — this is an indication that we should expect some finite-size 
dependence on our determination of the colloid-gas branch of the coexistence curves. The 
results reported in the next section confirm these expectations. 

Finally, we also report the distributions for the case / = 0.7, Rais/Rc = 1, '7p ~ 1-82 
[zpR^ = 0.85) (the average colloid-number histograms are reported in Fig. [2]), for which 
we have been unable to determine the coexistence fugacity using the equal-area method. 
The distribution of the colloid and polymer packing fractions at coexistence are reported in 
Fig. [5] and clearly explain the origin of the difficulties. The position of the colloid-gas branch 
varies significantly from sample to sample. We are thus far from the infinite- volume limit, 
since in this limit sample fluctuations are expected to disappear except close to the critical 
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FIG. 4: (Color online) Distribution of the coexisting colloid (top) and polymer (bottom) packing 
fractions for r]^ ^ 1.24 {zpR^ = 0.58), / = 0.4 and three different values of Rdis/Rc- 0.2, 0.6, 1.0. 
The colloid-gas phase data are reported in black, while the colloid-liquid phase data are reported 
in gray. Data from simulations with L/R^ = 16. 

point. These results provide again evidence that size effects are large for these values of the 
parameters. 

D. Finite-size effects 

The analyses presented above show that size effects may still be relevant for the data 
with L/Rc = 16. They appear to increase with increasing / and/or Rdis/Rc and should 
be particularly large for / = 0.7 and R^islRc = 1-0. To investigate size effects we have 
performed additional simulations with L/Rc = 20. In Fig. [6] we compare the results for the 
coexistence curve obtained by using these two different box sizes. We report results both in 
terms of rjc and rip and also in the reservoir representation (see insets) in terms of rjc and 1]^. 

On top we show the results for / = 0.4. Corrections here appear to be under control: the 
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FIG. 5: Distribution of the coexisting colloid (left) and polymer (right) packing fraction for rjp 
1.82 {zpRl = 0.85), / = 0.7 and i?dis/^c = 1- The colloid-gas phase data are reported in black, 
while the colloid-liquid phase data are reported in gray. Data from simulations with L/Rc = 16. 

plots in the terms of rjc and rjp show a small size dependence, while those in the reservoir 
representation appear to be reliable except close to the critical point, which is not unexpected 
since size corrections are large in a neighborhood of a second-order phase transition. For 
/ = 0.7, the colloid-liquid branch is determined quite reliably. On the other hand, the colloid- 
gas phase boundary varies significantly when L changes, especially for the case i?dis/ Rc = 1-0. 
This is not unexpected, given the results shown in the previous sections. Indeed, in all cases 
the polymer and colloid histograms are characterized by very narrow colloid-liquid peaks 
whose positions have a tiny dependence on the fugacity Zc, so that, even if z* is not precisely 
determined, the determined values rjcUq and r/p uq are quite reliable. In the colloid-gas phase 
the distributions are instead very broad, a clear indication that we are far from the infinite- 
volume limit. As we have already remarked, it is also possible that, for some values of the 
parameters, the system is in the one-phase region, even if the finite-size data show double 
peaks. 



15 



0.8 
0.6 
0.4 
0.2 
0, 



▲ 


1.2 


-▲ ▲ 


- 

▲ 
Aa 


1.1 
1 


A ^ 


▲ 


0.1 0.2 0.3 


A (0.4,0.20)L16 
w (0.4,0.20)L20 




▲ 







0.1 0.2 



0.3 




0-1 A (0.4,1.00)L16 
^ (0.4,1.00)L20 



0.05 0.1 0.15 0.2 0.25 



0.6 
0.5 
0.4 
0.3 
0.2 
0.1 






A 




A 


1.75 


A 
A 


■ 


A 
ft 




~A 




A 


^ 1.70 




■ 


■ 




- A 




A 


1.65 


- A 




A 






A 


A 


1.60 




1 A 1 


,A 



0.1 0.15 0.2 0.25 
lie 



A (0.7,0.20)L16 
■ (0.7,0.20)L20 



0.1 



0.15 



0.2 



0.25 



0.5^ 



0.4 -\ 



0.3- 

0.2- 

0.1^ 
▲ 



A 

A. 



1.8 
1.7 
1.6 
1.5 



k. U 



J , L 



0.1 0.15 0.2 0.25 



(0.7,1.00)L16 
(0.7,1.00)L20 

Ol5 



0.2 



0.25 



FIG. 6: Fluid-fluid binodal curves for / = 0.4, Rdis/Rc = 0.2 (top left), / = 0.4, Rdis/Rc = 1-0 (top 
right), / = 0.7, Rdis/Rc = 0.2 (bottom left), / = 0.7, Rdis/Rc = 1-0 (bottom right). We report the 
results for L = l6Rc and L = 20Rc in the 7]c, r]p plane and (inset) in the reservoir representation 



E. Demixing curves in the reservoir representation 

The estimates of [we report the average of z*{c) and zl{p)] as a function of r/^ for 
L/Rf. = 16 are reported in Fig. [71 In order to compare our results with those of Ref. 45| . 



we also report the estimates of the polymer reservoir packing fraction 77" at coexistence in 
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FIG. 7: Estimates of z*R^ as a function of the polymer reservoir packing fraction r]p (left) and of rjp* 
at coexistence in terms of the colloid reservoir packing fraction rj^ (right). Results for L/Rc = 16. 
In the legend the first number corresponds to /, while the second one gives the ratio Rdis/Rc- 



terms of the colloid reservoir packing fraction 77^ 58[. Note that, on a logarithmic scale, 
the values z*Rl for each / and Rdis/Rc lie quite precisely on a straight line, indicating 
that the colloid chemical potential at coexistence is well approximated by a linear function 
in rjp. Moreover, the position of the demixing curve depends essentially only on /. The 
ratio Rdis/Rc, hence the topological structure of the matrix, does not change significantly 
the coexistence curve. We have not performed a careful finite-size scaling analysis close to 
the critical point (a detailed discussion of the methods appropriate for random-field Ising 



critical points is reported in Refs. 



47 



59, 



601 ]) and thus we are not able to estimate //^ crit 
and 77c.crit precisely. We only note that for L/Rc = 16 and L/Rc = 20 double peaks are 
observed only for rj^ > 1.00, 1.03 for all three values of Rdis/Rc- We can thus set the lower 
bounds r^p crit ^ 1-03 and r^c.crit ^ 0.38. For / = 0.7 size effects are significantly larger than 
for / = 0.4, but we can still obtain the bounds ^^p ^rit ~ 1-6) '7c,crit ~ 0.405. 

We can use our results to understand qualitatively the behavior of a bulk colloid-polymer 
mixture in chemical equilibrium with the same dispersion adsorbed in a porous matrix. The 
main question is whether one can observe different phases in the bulk and in the matrix. If 
we use T]^ as control parameter, we see that for rj^ < "/^ccritbuik ~ 0-37 there is no transition, 
neither in the bulk nor in the matrix. If rj^ is larger, one may have a transition in the 
bulk and no transition in the matrix, given that ^^ccrit increases with /. For instance, for 
/ = 0.4 and 0.37 < r^^' < 0.385 we only observe phase separation in the bulk. If we increase 
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FIG. 8: Binodal curves (obtained by interpolating the data reported in Fig. [6]) for / = 0.40 (left) 
and / = 0.70 (right). In both plots we also report the bulk binodal curve and the effective critical 
points. In the legend the first number corresponds to /, while the second one gives the ratio 
Rdis/Rc- Simulations in a box of size L = WRc- 

further rj"^ [rjl^ > 0.385 for / = 0.4) the behavior is more complex, since phase separation 
occurs both in the bulk and in the matrix. For small rj^ the mixture is in the colloid-liquid 
phase both in the bulk and in the matrix. If rj^ is increased, that is polymers are added, 
the bulk coexistence curve is reached, see Fig. [71 Above the demixing line, one observes two 
different phases: in the bulk the system is in the colloid-gas phase, while in the matrix a 
colloid-liquid phase occurs. Thus, the presence of the matrix may induce, for certain values 
of the parameters, a capillary condensation of the colloids. Finally, for large rj^ above the 
matrix coexistence curve, a colloid-gas phase occurs both in the bulk and in the matrix. 



F. Binodals in the system representation 

In Fig. |8]we report the results for the binodals (we interpolate the MC data for L/Rc = 16) 
as a function of rjc and rjp. For / = 0.4 they are quite close to the bulk binodal curve and show 
only a tiny dependence on the ratio Rdis/Rc- In the figure we also report an estimate of the 
critical point obtained by determining the intersection of the diameter with the interpolation 
of the coexistence data. This provides a very rough estimate of the critical parameters, which 
can only be accurately determined by performing a careful finite-size scaling analysis. In 
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the bulk the analysis of the results with L/Rc = 16 givesjjccrit ~ 0.12, ?7p,crit ~ 0.31, which 
should be compared with the precise determination 3l|, |32| 



r/e,crit = 0.1340(2), 77p,crit = 0.3562(6). (17) 

Apparently our simple extrapolation underestimates ?7c,crit by 10% and r^p^crit by 15%, which 
we can take as indications of the systematic error. If we perform the same analysis for / = 0.4 
we obtain recent ~ 0.17 for all values of Rais/Rc'- the colloid packing fraction at criticality 
is essentially independent of the matrix topology and increases with increasing /. As for 
the polymer packing fraction, the dependence on the size ratio Rais/ Rc is somewhat larger. 
The value r^p^crit decreases with increasing R^is/Rc and varies between 0.23 and 0.29. If we 
assume that the results for L/Rc = 16 underestimate the correct, infinite- volume results by 
10% and 15% as they do in the bulk, we would guess ?7c,crit ~ 0.19 and 0.27 < r^p^crit ^ 0.34. 

As we have already stated our results for / = 0.7 are only reliable in the colloid-liquid 
phase. In this regime the binodal curve depends somewhat on Rdis/Rc- for Rdis/Rc = 0.2 it 
is close to the bulk curve, while for R^is/Rc = 1 it is significantly below it. The critical point 
position is consistent with what was observed for / = 0.4: recent shows little dependence on 
Rdis/Rc, while '7p,crit decreases with increasing R^is/Rc- 



We can compare our results with those obtained in Refs. [45|, |48|. By using density- 



functional methods, Ref. [45|] studied the model with Rdis/Rc = 1 at the slightly lower value 
of /, / = 0.37 (corresponding to r/dis = 0.05), and several values of g, g = 0.3, 0.6, 1.0 
(none of them agrees unfortunately with ours). Their results are not consistent with ours. 
They find that in all cases the binodal curve in the presence of the matrix is below that in 
the bulk, while here we find the opposite except deep in the colloid-liquid region, i.e. for 
rjc > 0.20. Moreover, they find that ?7c,crit decreases with increasing / in all cases, which 
is again in contrast with our results. On the other hand we are fully consistent with the 
results of Ref. [48|, which study the model for Rdis/Rc = 1, / = 0.37 (?7dis = 0.05), and 
q = 0.8. For the critical point they obtain ?7c,crit = 0.192 and r^p^crit = 0.292, confirming 
that ?7c,crit increases in the presence of disorder. From a quantitative point of view, their 
critical-point estimates are fully consistent with ours. In particular, the naive extrapolation 
we have performed above apparently estimates correctly the critical-point position at the 
5% level. 
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IV. CONCLUSIONS 



In this paper we determine the fluid-fluid demixing curves for the AO model in a colloidal 
matrix for q — 0.8 and several values of Rdis/ Rc and / (equivalently, of the disorder concen- 
tration Cdis)- This study should provide quantitative informations on the phase behavior of 
a polymer-colloid mixture in a porous material close to the 9 point. Our main results are 
the following: 

• Disorder is specifled by two parameters: the disorder packing fraction r/dis and the 
ratio Rdis/Rc- At least for R^is/Rc < 1, the parameter range we consider, most of the 
disorder dependence of the results can be parametrized by using a single parameter, 
the effective matrix fiUed-space ratio /, which gives the volume fraction unavailable 
to the colloids. In the Zp, Zc plane (or equivalently in terms of the reservoir packing 
fractions 7)1 and r^p the coexistence curve depends essentially only on /. 

• It is possible to observe capillary condensation of the colloids. For certain values of 
the parameters a colloid-gas bulk phase is in equilibrium with a colloid-liquid phase 
in the matrix. 

• At least for / < 0.4 the binodal curves expressed in terms of the packing fractions 
(system representation) show a relatively small dependence on disorder. The critical 
point instead changes signiflcantly. The critical colloid packing fraction 77c,crit is, to a 
large extent, only a function of / and it increases as / increases. The critical polymer 
packing fraction rjp^crit depends instead both on / and Rdis/Rc- At fixed / it decreases 
as Rdis/Rc increases. 

The authors gratefully acknowledge extensive discussions with Ettore Vicari. The MC 
simulations were performed at the INFN Pisa GRID DATA center and on the INFN cluster 
CSN4. 

Appendix A: Monte Carlo simulations: some technical details 

We have performed simulations in the grand-canonical (GC) ensemble, which physically 
describes a system adsorbed in a colloid matrix in chemical equilibrium with a reservoir of 



20 



pure polymers and a reservoir of pure noninteracting colloids. The basic parameters are the 
colloid and polymer fugacities and z^. In the bulk the partition function is 

z,, z,) = J2 ^p'^c'QiV, N„ AT,), (Al) 

where QiV, Np, Nc) is the configurational partition function of a system of Np polymers and 
Nc colloids in a volume V. We drop the irrelevant thermal length so that Q{V,1,0) = 
Q{V, 0, 1) = V. In the presence of first-order transitions it is quite difficult to sample 
correctly the GC distribution. To bypass the difficulties we use the umbrella-sampling 
(sometimes also called multicanonical) method jssl. Instead of generating configurations 
with the GC weight, we use an umbrella distribution 

^ z^^z^^e-^^ (A2) 

with a properly chosen tt{Nc) which is defined below. If {■)gc and (■),r are the averages with 
respect to the GC distribution and to the distribution (1A2I) . respectively, we have 

(0(iV..iV,))„.^ <"^-";(,^f>»- . ,A3) 

This relation allows us to obtain GC averages from simulations using the distribution ( ]A2p . 
The function tt{Nc) must be chosen so that in the simulation the system can move easily 
between the two phases. Consider the histogram of Nc in the GC distribution, i.e. 

h{Nc,o) = {6{Nc,Nc,o))GC, (A4) 

where S{x,y) is Kronecker's delta. Assume that the system is close to phase separation so 
that h{Nc) has two peaks at A^c,min (colloid-gas phase) and A^^max (colloid- liquid phase). The 
optimal choice is then 

7l{N) = ah{Nc,min) N < Nc,mi^, 

7c{N) = ah{N) Nc,min <N< AT,,^^, (A5) 

n{N) = ahiNc,^^,) N > 



' c,max5 



where a is an irrelevant constant. Indeed, if A"c,min ^ < A^c,max; the observed histogram in 
the umbrella distribution is fiat, i.e. independent of N^. Hence, the system can move freely 
between the two phases, allowing a precise determination of any required thermodynamic 
property. 
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In order to determine the colloid fugacity z* at coexistence for a given value of the 
polymer fugacity Zp, we consider Nm colloid fugacities {zc,m}, such that for 2:^1 (-2c,Af™) the 
system is in the colloid-gas (colloid-liquid) phase. Then, we determine the umbrella functions 
T^miNc) iteratively. First, we perform a short hysteresis cycle in which we perform A^therm 
GC iterations at Zc = -Zc.i, then at Zc^2, and so on, up to Zc^n^; then we decrease Zc till we 
reach again 2:^1. If hm'~^{Nc) and hm' (Nc) are the histograms obtained at z = Zm (the 
+ refers to the distribution obtained while increasing Zc and the — to that obtained while 
decreasing the fugacity), we set hm{Nc) = hm'~^{Nc) + hm~{Nc) and 

7r£^ ( AT,) = hil^ {N,)/M if /i^) ( AT,) > M, 

(A6) 

where M = ma.XN^[hm\Nc)]/10. Then, we repeat again the same hysteresis cycle several 
times. At iteration k, for each Zc,m we perform the simulation using the distribution (IA2p 
with TT = 7Cm ■ Then, we set 

7Ti!:\N,) = n^J:''\N,)hl^\N,)/M if hl^\N,) > M, 
n^J:\N,) = 7ri!:-'\N,) if hi^\N,) < M, 

where M = maxN^\hm\Nc)]/10. We stop when we observe that, for at least some values 
of m, hin^'^iNc) and hl^^'~{N,) are nonvanishing in an interval of values of Nc that extends 
between the two phases. 

Once we have a reasonable estimate of the functions UmiNc), we could just perform an 
extensive simulation at single value of Zc^m, (an optimal choice would be to take the value for 
which TCmlNc) is clearly bimodal). Data for different values of Zc could just be obtained by 
standard reweighting techniques. However, we have found more convenient, to use all infor- 
mation we have collected and simulate all systems together, using the simulated-tempering 
method (sgI. Note that, in the standard implementation of the method, one should be care- 
ful that the fugacities Zc^m are such that the colloid-number distributions overlap; otherwise, 
no fugacity swap is accepted. In our case, since we use umbrella distributions, the overlap 
condition is always verified, and thus the number A^ of needed systems is always small. 
Typically we take Nm = 10. If 

^n^{V, Zp, Zc,m) = Yl ^nfr^^^' ^^)' (^^) 

22 



we consider the extended partition function 



The constants fm are chosen so that all terms in the sum are approximately equal. If we 
require 



we obtain 



with 



Rm ^ = S^, (All) 



m—1 Jn 



Rm ^ r-^) , (A12) 















^€,171—1 / 


1 ttUNc) 



TT,m 



7r,m— 1 



where (■)7r,m indicates the mean value with respect to the umbrella distribution ( IA2p with 
Zc = Zc,m, TT = TTm- Combining these expressions we define the ratios as 

= ^/RJs;::. {au) 

Jm— 1 

The constants Rm and Sm are determined together with the umbrella sampling functions 
7im- Then, we set /i = 1 and use Eq. flA14p to determine the constants fm, ^ > 2. 

In the matrix case, the GC partition function is still given be Eq. f lAip . with the only 
difference that one should take into account the interactions between the freely moving 
particles and the matrix. Since the GC partition function depends on the matrix, also the 
functions vr^ and the constants fm are matrix dependent. Thus, we recompute them when 
we restart the simulation with the different matrix. 

In the MC simulations we take Nm ~ 10. One MC iteration consists in 3 fugacity swaps 
and 1000-5000 GC moves in which colloids and polymers are inserted or removed. For this 



purpose we use the cluster move discussed in Ref. [57| together with standard moves in 
which a single polymer is removed or inserted. For each disorder instance, we perform A^inj 
iterations to determine the umbrella functions and then A^iter iterations to measure several 
histograms. Typically, Aini varies between SOOOA'^ and 20000A^m, while Alter is of the order 
of 20000Ar^. 
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In the simulation we determine the colloid and polymer histograms for a large number 
(typically 100) of colloid fugacities Zc,r- They are obtained by measuring, for each matrix 
realization a, the reweighted histograms Pc(«, -2c,r, Nc,o) and Pp{a, Zc,r, Np,o)- 

Pci<y,Zc,r,Zc,m,Ncfl) = ^1 ( ) ^m{Nc,i)5iNc,i, Ncfl)5 iZc,m, Zc,i) , (A15) 

\ Zc m I 

Pp{a,Zc,r,Zc,m,Npfl) = X] ( ) TT^ ( A^c,i)5 (iVp.i , A/p,o)5 (^^c,™, ^c,j) , (Al6) 



where i refers to the MC iteration, and Np^i, N^^i, Zc^i are the number of polymers and colloids 
and the colloid fugacity at the ith iteration. The colloid histogram is then 



hc,ave{Nc 0, Zp, Zc r) ~ nr / , ST^ / ]\t \ 

Na ^ I2_^j^ Pc{a, Zc,r, Zc,m, Nc) 



(A17) 



where is the number of matrix realizations. Note that we obtain a different estimate of 
the distributions at Zc^r for each of the Zc,m- 

As a final comment, note that our estimates ( 1A17P are biased, since they are disorder 
averages of a ratio of thermal averages. This means that, if we take the limit Na — )■ oo 
at fixed A^iter, we obtain estimates that differ from the correct result by a term (the bias) 
of order 1/A^itcr- One could perform a bias correction, as discussed in Ref. [6l|. However, 
given the small number of disorder instances, we have found that in the present case the 
bias correction is not relevant. 
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